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Abstract 

A piecewise-deterministic Markov process is a stochastic process whose behavior is governed by an 
ordinary differential equation punctuated by random jumps occurring at random times. We focus on 
the nonparametric estimation problem of the jump rate for such a stochastic model observed within a 
long time interval under an ergodicity condition. We introduce an uncountable class (indexed by the 
deterministic flow) of recursive kernel estimates of the jump rate and we establish their strong pointwise 
consistency as well as their asymptotic normality. We propose to choose among this class the estimator 
with the minimal variance, which is unfortunately unknown and thus remains to be estimated. We also 
discuss the choice of the bandwidth parameters by cross-validation methods. 
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1 Introduction 


Piecewise-deterministic Markov processes (PDMP’s in abbreviated form) have been introduced in the 
literature by Davis in [17] as a general class of continuous-time non-diffusion stochastic models, suitable 
for modeling deterministic phenomena in which the randomness appears as point events. The motion of 
a PDMP may be defined from three local characteristics: the flow the jump rate A(x) and the 

transition measure Q{x,Ay). Starting from some initial value Xq, the process evolves in a determinis¬ 
tic way following <l?(Xo,t) until the first jump time Ti which occurs either when the flow reaches the 
boundary of the state space or before, in a Poisson-like fashion with non homogenous rate A(<l?(Xo, t)). 
In both cases, the post-jump location of the process at time Ti is governed by the transition distribution 
Q{^{Xo, 21), dy) and the motion restarts from this new point as before. This family of stochastic models 
is well-adapted for tackling various problems arising for example in biology [9, 16, 30, 31, 32, 33, 34], 
in neuroscience [22] or in reliability [7, 14, 13, 18]. Indeed, most of applications involving both de¬ 
terministic motion and punctual random events may be modeled by a PDMP. Typical examples are 
growth-fragmentation models composed of deterministic growths followed by random losses. For exam¬ 
ple [33] , the size of a cell grows exponentially in time, next the cell divides into two offsprings whose size 
is about half size of the parent cell, and so on. 

Proposing efficient statistical methods for this class of stochastic models is therefore of a great interest. 
Nevertheless, the very particular framework involving both deterministic motion and punctual random 
jumps imposes to consider specific methods. For instance, the authors of [5] have shown that the well- 
known multiplicative intensity model developed by Aalen [1] for estimating the jump rate function does 
not directly apply to PDMP’s. Alternative approaches should be thus proposed. In the present paper, 
we focus on the recursive nonparametric estimation of the jump rate of a PDMP from the observation of 
only one trajectory within a long time interval. More precisely, the purpose of this work is to show how 
one may obtain by kernel methods a class of consistent estimators for the jump rate, and how one may 
choose among this class in an optimal way. To the best of our knowledge, the nonparametric estimation 
of the jump rate in a general framework has never been investigated. 

As PDMP’s may model a large variety of problems, some methods have been developed by many 
authors for their statistical inference. As presented before, the randomness of a PDMP is governed by 
two characteristics: the transition kernel Q{x,dy) and the jump rate A(x). As a consequence, two main 
questions arise in the estimation problem for such a process, namely the statistical inference for both 
these features. On the one hand, a few papers investigate some nonparametric methods for estimating 
the transition function of a PDMP either for a specific model [13] or in a more general setting for a 
d-dimensional process [3]. On the other hand the estimation of the jump rate A( 2 ;) or of the associated 
density function has been more extensively studied by several authors. Without attempting to give an 
exhaustive survey of the literature on this topic, one may refer the reader to [5, 19, 20, 25, 28] and the 
references therein. In the book [25], the author studies likelihood processes for observation of PDMP’s 
which could lead to inference methods in a parametric or semi-parametric setting. The papers [19, 20] deal 
with the nonparametric estimation for some PDMP’s used in the modeling of a size-structured population 
observed along a lineage tree. In both these articles, the authors rely on the specific form of the features 
of the process of interest in order to derive the asymptotic behavior of their estimation procedure. These 
techniques have been generalized in [28] to introduce a nonparametric method for estimating the jump 
rate in a specific class of one-dimensional PDMP’s with monotonic motion and deterministic breaks, 
that is to say when the transition measure Q{x,dy) is a Dirac mass at some location depending on x. 
The procedures developed in these papers [19, 20, 28] are obviously of a great interest but strongly use 
the particular framework involved in the investigated models and are thus not well adapted in a more 
general setting. In [5], the authors show that the famous multiplicative intensity model only applies 
for estimating the jump rate of a modified version of the underlying PDMP. This leads to a statistical 
method for approximating the conditional density associated with the jump rate for a process defined on 
a bounded metric state space. 

A main difficulty throughout the present paper and the articles [3, 5] is related to the presence of 
deterministic jumps when the path reaches the boundary of the state space. This feature is often used for 
modeling a deterministic switching when the quantitative variable rises over a certain threshold [2, 24]. 
In a statistical point of view, the interarrival times are therefore right-censored by a deterministic clock 
depending on the state space, which leads to some technical difficulties. We would like to emphasize that 
the techniques developed in the references [19, 20, 25, 28] do not take into account the likely presence of 
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forced jumps in the dynamic. 

One may also find in the literature a few papers [6, 8, 10, 11] which focus on the estimation of various 
functionals for this family of stochastic models. More precisely, the authors of [10, 11] provide numerical 
methods for the expectations and for the exit times of PDMP’s. In addition, the article [8] deals with a 
PDMP introduced for modeling the temporal evolution of exposure to a food contaminant and consider 
a simulation-based statistical inference procedure for estimating some functionals such as first passage 
times. Plug-in methods have been studied in ]6] for a non-ergodic growth-fragmentation model which is 
absorbed under a certain threshold. In many aspects our approach and these papers are different and 
complementary. Indeed they are devoted to the estimation of some functionals of PDMP’s while we focus 
on the direct estimation of the primitive data of such a process. 

In this article, we introduce a three-dimensional kernel estimator computed from the observation 
of the embedded Markov chain of a PDMP composed of the post-jump locations and the travel 
times Sn+i along the path 3>(Z„,t). We establish its pointwise consistency as well as its asymptotic 
normality in Theorem 3.3. The estimate that we consider is recursive: it may be computed in real-time 
from sequential data, which may be relevant in many applications. We deduce from this result two 
first corollaries about the nonparametric estimation of the conditional density f{x, t) of the interarrival 
time Sn+i at time t conditionally on the event {Zn = x} (see Corollary 3.8) and its survival function 
G{x,t) = P(S'„+i > t \ Zn = x) (see Corollary 3.9). We also investigate in Corollary 3.10 the asymptotic 
behavior of an estimator for the composed function Ao ^{x,t) obtained as the ratio f{x,t)/G{x,t)- We 
derive in (18) an uncountable class (indexed by the states ^ hitted by the reverse flow 'I>(a;, —t) for some 
time t) of consistent estimates of the jump rate \{x). In other words, for each ^ = 'l>(a;, —t), we get 
a good estimate of \{x). We show how one may choose among this class of estimators by minimizing 
their asymptotic variance. We state in (19) that this procedure is equivalent to maximize the criterion 
along the curve ^{x,—t), i.e. ^ = ^{x, where iZoo(C) denotes the invariant 

measure of the post-jump locations Zn and Tx{^) is the only deterministic time to reach x following 
$(^,1). The choice of this criterion is far to be obvious without precisely computing the limit variance 
in the central limit theorem presented in Corollary 3.10. Indeed, a naive criterion to maximize is the 
invariant distribution Voa{^) along 4>(a;, —t): the larger i^oo(C) is, the larger the number of data around 
^ is and the higher the quality of the estimation should be. Nonetheless, this simple criterion does not 
take into account that the estimate also depends on the time of interest Tx{^) (see Remark 3.11). This 
question is also investigated from a numerical point of view in Subsection 4.2 where we show on synthetic 
data that the choice of the criterion Ti;(0) *s better than the naive one UoaiO- The bandwidths 

in kernel methods are free parameters that exhibit a strong influence on the quality of the estimation. 
We discuss the choice of the bandwidth parameters by a classic procedure that consists in minimizing the 
Integrated Square Error, computed here along the reverse flow 4>(a:, —t): we introduce a cross-validation 
procedure in this Markov setting and we prove its convergence in Propositions 3.14 and 3.15. Finally, we 
would like to highlight that the regularity conditions that we impose are non restrictive. In particular, 
neither the deterministic exit time from the state space is assumed to be a bounded function, nor the 
transition kernel is supposed to be lower-bounded, as is the case in [5] (see eq. (2) and Assumptions 2.4). 
In addition, the forms of the transition measure and of the deterministic flow are not specified. 

The sequel of the paper is organized as follows. We begin in Section 2 with the precise formulation of 
our framework (see Subsection 2.1) and the main assumptions that we need in this article (see Subsection 
2.2). Section 3 is devoted to the presentation of the statistical procedure and the related results of con¬ 
vergence. More precisely, a three-dimensional kernel estimator for the inter-jumping times is introduced 
and investigated in Subsection 3.1. We derive a class of estimators of the jump rate and propose how 
to choose among it in Subsection 3.2. The crucial choice of the bandwidth parameters is studied in 
Subsection 3.3. Finally, a self-contained presentation of the whole estimation procedure is provided in 
Subsection 4.1, then illustrated in the sequel of Section 4 on three different application scenarios, with 
various sample sizes and state space dimensions, involving both simulated and real datasets. More pre¬ 
cisely, we focus on the TCP window size process used for modeling data transmission over the Internet in 
Subsection 4.2. Estimation of bacterial motility is tackled in Subsection 4.3, while acceleration of fatigue 
crack propagation is considered in Subsection 4.4. The proofs and the technicalities are postponed in 
Appendix A, B and C at the end of the paper. 


3 


2 Problem formulation 


This section is devoted to the definition of the PDMP’s and to the presentation of the main assumptions 
that we impose in the paper. 

List of notations: In this paper denotes the Borel algebra of R'* endowed with the Euclidean 

norm | • |. In addition, the Lebesgue measure on S(R'*) is denoted by Ad(da;), with the particular notation 
in the one-dimensional case Ai(dx) = d®. The ball of R'* with radius r and center x is denoted by 
Bd{x, r). The closure of a set E is denoted by E while dE stands for its boundary. 

2.1 Definition and notation 

The motion of a PDMP may be described as the solution of an ordinary differential equation $ punctuated 
at random times by random jumps governed by a transition kernel Q (see Figure 1). The random jumps 
occur either when the deterministic motion hits the boundary of its state space, or before, with non 
homogeneous rate A taken along the curve defined by the differential equation <1?. 



Figure 1: Schematic path of a PDMP. This paper is devoted to the nonparametric estimation of the state- 
dependent rate A, that governs the spontaneous generation of jumps, from a long-time trajectory. 

More precisely, a PDMP (Xt)t>o on is defined from its three local characteristics (A, 

• 3> : X R —>■ is the deterministic flow. It satisfies the semigroup property, 

V® € R"^, Vt, s € R, 'l>(®,t -I- s) = t), s). 

• A : R'* —>■ R+ is the jump rate. 

• Q : (R'*, f3(R‘*)) —>■ [0,1] is the transition kernel. 

We define the deterministic exit times of E for the flow and for the reverse flow as, for any x £ E, 

= inf {t > 0 : ^{x,t) £ dE} and (®) = inf {t > 0 : ^{x,—t)£dE}. (1) 

In all the sequel, we consider a PDMP (Xt)t>o evolving on an open subset E of R'*. In this context, 
we impose as usual [17, (24.8) Standard conditions] that, 

3e>0, f A(<l>(®, t)) dt <-|-(XJ, 

Jo 

and 

yx£E,Q{x,E\{x}) = l. (2) 

In addition, we restrict ourselves to the case where the transition kernel Q admits a density with respect 
to the Lebesgue measure, 

yx£R‘^, yA£B{R‘^), Q{x,A)= f Q{x,y)Mdy). (3) 

J A 

This assumption is natural when one considers multivariate real-valued PDMP’s, and is satisfied in 
various problems arising in biology [9], population dynamics [6] or in insurance [17, (21.11) An insurance 
model]. 
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Starting from any initial condition Xq = x, the motion of {Xt)t>o may be described as follows. The 
distribution of the first jump time Ti is given by, 

Vf > 0, P(Ti >t\Xo = x) = h^^{- fo «)) ds) if f (4) 

0 else. 

In other words, the process jumps either when the flow hits the boundary of the state space at time t~^(x) 
or in a Poisson-like fashion with rate A o $ before. Next the post-jump location Zi at time Ti is defined 
through the transition kernel Q: for any test function ip, we have 

E lip(Zi) I Ti, Xo = x] = J ip{u) Q{^{x,Tt_),Au). 

The path between 0 and the first jump time Ti is given by, 


VO < t < Ti, Xt = 


cf>(x, t) if t < Ti, 
Z\ else. 


Now starting from the post-jump location Xt^ , one chooses the next inter-jumping time S 2 = T 2 — T\ 
and the future post-jump location Z 2 in a similar way as before, and so on. One obtains a strong Markov 
process with {Tn)„>o as the sequence of the jump times (where To = 0 by convention). The inter-jumping 
times are defined by So = 0 and, for any integer n > 1, Sn = Tn — r„_i. Finally (Z„)„>o denotes the 
stochastic sequence of the post-jump locations of (Xt)t>o, with for any n, = Xt„. 

All the randomness of the PDMP (Xt)t>o is contained in the stochastic sequence (Zn, S„+i)n>o 
which is a Markov chain. In addition, the post-jump locations (Zn)n>o also form a discrete-time Markov 
process on the state space E because of the condition (2). In this paper, denotes the distribution of 
the n**' post-jump location for any integer n > 0, while V denotes its Markov kernel, 

Vx€E, VA€l3(R‘'), P(x,A) = P(Z„+i e A\Z„ = x) 

= f S{x,dt) Q{^{x,t),A), (5) 

ilR^. 


where S stands for the conditional distribution of S„+i given Zn. for any n, 

M X £ E, Vt > 0, 5(x, (t, -foo)) = P(S'„+i > t I Zn = x) 

^ I exp(^-/Q A('F(x,s))ds) At<t+{x), 

( 0 else, 

in light of (4). We would like to highlight that the conditional distribution Six, ■) is absolutely continuous 
with respect to the unidimensional Lebesgue measure on (0, t^(x)) with sometimes a singular component 
at f^ix), 

VAg5(IR+), 5(x,A)= [ f{x,t)dt + 5(x,An{t+(x)}), (7) 

J An(0,t+(a;)) 

where the conditional density / may be obtained by deriving (6), 


y X £ E, y 0 < t < f^ix), fix, t) = A(tE>(x, t)) exp ( — / A(4?(x, s))ds ) . 


(-r 


In all the sequel G stands for the conditional survival function associated with /, that is, 
yx£E, yQ<t<t'^ix), G(x, t) = 5(x, (t,-boo)). 


( 8 ) 

(9) 


where S is defined in (6). As highlighted before, the process (Zn, S'n+ijn^o forms a Markov chain on the 
set F defined by 

T= (J {x}x [0,t+(x)]. (10) 

x£E 

TZ denotes the transition kernel of this process. 


V(x,t) £F,yAxB£ Z3(R‘' x R+), 7^((x,^),A x B) 


P(Zn-t-l £■ A, Sn+2 £■ B I Zn — X, Sn + 1 — t) 

f Q(<&(x,t),d05(C.B), (11) 


and, for any n, pn denotes the distribution of the couple (Z„, 5'„+i). 
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2.2 Assumptions 

The main assumption that we impose in the present paper is a condition of ergodicity on the Markov 
chain {Zn)n>o- This property is often a keystone in statistical inference for Markov processes and may 
be directly imposed [3, 5] or established [28] from the primitive features of the data. 

Assumption 2.1 There exists a distribution Vao on E such that, for any initial distribution vq = 

X £ E, 

lim ||l/„ — Vao'^TV = 0, 

n—^ + oo 

where || ■ ||tv stands for the total variation norm. 

This assumption may be checked directly on the Markov kernel V of {Zn)n>o from the existence of a 
Foster-Lyapunov’s function or Doeblin’s condition for instance [29, Theorem 16.0.2]. In the following 
remark we establish a first property of the sequence (i/'n)n>o and of its limit Uoa- 

Remark 2.2 Since the transition kernel Q is assumed to be absolutely continuous with respect to the 
Lebesgue measure (3), the kernel V given by (5) of the post-jump locations {Zn)n.>o olso admits a density. 
As a conseguence, for any integer n, the distribution of Z^ and thus the invariant measure 
introduced in Assumption 2.1 admit a density on the state space E. For the sake of clarity, we write 
Vooidx) = i'oa{x)\d{dx) with a slight abuse of notation. 

We add some regularity conditions on the main features of the process to show the convergence of the 
estimates in Theorem 3.3. 

Assumptions 2.3 

• The sup-norms |jQ||oo and ||/||tx> are finite. These conditions are used in the proof of Theorem 3.3 
to find an upper bound of the non diagonal terms of some square variation process of interest. 

• The functions Q and f are Lipschitz, 

yx,y £ E, y0< s <t< t+(a;) At+(y), \f{x,t) - f{y,s)\ < [f]Lip{\t - sj + \y-x\), 

'ix £ E, \/y, z £ E,\Q[x,y) - Q{x,z)\ < [Q]up\y - z\. 

These conditions are used in the proof of Theorem 3.3 to control the diagonal terms of the same 
variation process and to study the convergence of some remainder terms. 

• The survival function G is Lipschitz, 

Vx, y £ E, VO < t < t'^(x) A IG(x,t) - G(y,t)l < [G]Lip\x - y\. 

This condition is used in the proof of Theorem 3.3 to investigate the convergence of the remainder 
terms. 

• The deterministic exit time is continuous. This condition is used to find some admissible initial 
bandwidths vq and Wq. 

Finally, we consider an additional condition on both the transition kernels V and Q and the flow $ in 
order to ensure the Lipschitz mixing property of the Markov chain {Zn)n>o- This will be sufficient to 
establish the almost sure convergence to 0 of the remainder term with the adequate rate in the proof of 
Theorem 3.3. 

Assumptions 2.4 The transition kernel V of the Markov chain (Zn)n>o satisfies, for some Ui > 1 and 

02 < 1 , 

V{x,y)£E‘^, f \u — v\°'^V{x,du)V{y,dv) < a 2 \x — y\°‘^. 

In addition, the composed function Q(4>(-,'),-) belongs to the regularity class Li(ri,r 2 ) defined by 
\Q{^{xi,ti),yx)-Q{tI>{x 2 ,t 2 ),y 2 )\ = O {\{xi,ti,yi) - {x 2 ,t 2 ,y 2 )P {\{xi,G,yi)\"^ -I- \{x 2 ,t 2 ,y 2 )\"^ + 1)), 
for some positive numbers ri and r 2 satisfying 2(ri -|- r 2 ) < ai. 
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3 Estimation procedure 

3.1 Inference for the inter-jumping times 

For any integer n, we introduce the <7{Zo, Si, ..., Zn-i, S'n)-nieasurable functions 9^ defined on E, and 
JF" and defined on the interior F (let us recall that the set F is given in (10)), by, 


V {x, t) e F, F"{x, t) 

V {x,t) £ F, Q"{x,t) 

'ix £ E, 9^{x) 



( 12 ) 

(13) 


where Kp denotes a kernel function on R^, p € {l,d}, and the band widths are defined for any integer k 
by Vk = vo{k + 1)““ and Wk = wo{k + 1)“^ for some a, j3 > Q and initial positive values vo and wo- 
It should be already noted that these quantities are of a great interest in the statistical study of the 
inter-jumping times of the PDMP {Xt)t>o- Indeed, we will see that: 

• The ratio estimates the conditional density f{x,t) defined in (8). 

• The ratio estimates the conditional survival function G{x,t) defined in (9). 

• The ratio estimates the composed function X(^(x,t)). 

In addition, as the name suggests, we will state that 9^{x) is a good estimate of the density Vao{x) of 
the unique invariant distribution of the post-jump locations, which is relevant in the estimation problem 
for PDMP’s but has already been investigated in [3, Proposition A.11]. 

In all the sequel, we impose a few assumptions on both the kernel functions Ki and IKtj. 


Assumptions 3.1 For any p £ {l,d}, the kernel function Kp is assumed to be a nonnegative smooth 
function satisfying the following conditions: 

• The sup-norm ||Ap||oo is finite. 

• /rp Il'pdAp = 1. 

• suppKp C Bp{0p,S). Together with well-chosen initial bandwidths, this condition avoids to compute 
the kernel estimator from data located at the boundary of the state space F (see also Remark 3. 4). 

• Only for p = d. The function is Lipschitz, 

yX, y £ R'*, |Kd(a;) - IKd(j/)| < [Rjuipja; - y\. 

This condition is used to show eq. (57) in the proof of Proposition 3.14- 

Remark 3.2 In particular, Assumptions 3.1 ensure that, for any p £ {l,ci}, Tp = KpdAp is finite. 
This is used to find an upper bound of the non diagonal terms of the variation process in the proof of 
Theorem 3.3. In addition, the integral |u|Kp(u)dAp(u) is also finite, which is needed to establish the 
almost sure convergence of the remainder terms in the proof of the same result. 


In the sequel the admissible set for the bandwidth parameters a and (3 is given by 

A = {(a, /3) G : a > Q, fi > Q, ad -\- P < 1, ad -\- fl -\- 2min(a, (1) > 1} . 

In this part, our main result is obtained from the use of vector martingales and is stated in the following 
theorem. 


7 










Theorem 3.3 For any couple (x,t) € F such that Uao(x)f(x,t) > 0, for any (a,f5) € A and (wo,wo) 
such that 

t + wo5 < inf (14) 

where 5 appears in the third item of Assumptions 3.1, we have the almost sure convergence, 


■ F"(x,t) ' 


12 oo{x)f{x,t) 

g”(x,t) 

a.s^ 

Uoo{x)G{x,t) 



Voa (x) 


and the asymptotic normality. 


l-ad-/3 1 

■ T^{x,t) ' 


f{x,t) 

n 2 

g^{x,t) 

— 

Voo{x)G{x,t) 

V 

l J 


IXooix) 


d 


Af{03,T,{x,t,a,l3)), 


where the variance-covariance matrix T,{x,t, a, fd) is degenerate with only one positive term at position 
(1,1). T,{x,t,a,fd) is defined in (52). 

Proof. The proof is stated in Appendix B. □ 

Remark 3.4 The existence of a couple {vq,wo) satisfying (14) is obvious whenever the exit time is 
continuous (see Assumptions 2.3). This condition ensures that all the inter-jumping times used in the 
calculus of T'^ and are not obtained from forced jumps when the process reaches the boundary of the 
state space. In the case where T^ix) = oo, it is obvious that no interarrival times are right-censored. The 
consistency and the asymptotic normality are therefore still accurate without any condition on (vo,wo). 


Remark 3.5 In Theorem 3.3, the choice (14) of the initial bandwidths vo and wo is locally dependent 
on the point of interest. This may appear restrictive but may be avoided by considering the elements of 

C = |c X [0,r] : C is a compact subset of E and T < in^t'''(a;)j-. (15) 


Indeed, for any C x [0,r] G C, there always exists a couple (vo,wo) such that 

T-\-wqS < inf inf t^(£). 


Thus, (vo,wo) satisfies (14) for any point {x,t) £ C x [0, T]. 


Remark 3.6 The variance-covariance matrix appearing in the asymptotic normality presented in The¬ 
orem 3.3 is degenerate with only the component (1,1) positive. It means that the rate of the estimators 
vO, and Q" is faster than the one of F". This is straightforward because is obtained by smoothing 
the empirical distribution of the data both in the spatial and temporal directions contrary to and Q". 
The proof of the previous result may be adapted to show the two-dimensional central limit theorem, with 
a such that ad < 1 and a{d + 2) > 1/ 


l-ad 

n 2 


( 

■ g"ix,t) ■ 



Z2S,(x) 



Vac{x)G{x,t) 

Voo{x) 


N{Q 2 ,'E' {x,t,a)), 


where S'(x,t,a) is a diagonal 2 x 2-matrix. The keystone to state this convergence is the behavior given 
in (50) of the hook of a vector martingale. 


If we assume the geometric ergodicity of the Markov chain {Z„)n>Q, one may also obtain the rate of 
convergence of the variances of the estimates and uniformly on any compact subset of F such 

that the parameters vo and wo may be uniformly chosen (see Remark 3.5). 















Proposition 3.7 Let us assume that there exists 6 > 1 such that \\vn. — Voo\\TV = 0{h~"). The geometric 
ergodicity is in particular ensured by Doeblin’s condition (see [29, Theorem 16.0.2]). Then, for any set 
C X [0, T] G C and for any couple [a, (3) such that 2(ad + /3) < 1, we have 


sup Var(J^"'(a:, t)) 
(x,t)ec 


sup 

(x,t)€C 


(Va.r{g'^{x,t)) + Var(PSj(a:))) 


/ 2ad — 1 \ 

0(n ). 


Let us recall that C has been defined in (15). The rate of convergence for T'^ is faster than the one given 
in Theorem 3.3 whenever 3{ad + 13) < 1. 

Proof. The proof is similar to the demonstrations of Proposition B.5 and Corollary B.6 of [3] and relies 
on the control of the covariance process of functionals of a geometrically ergodic Markov chain (see [29, 
Theorem 16.1.5]). □ 

We present in the sequel some corollaries of Theorem 3.3 that are of interest in the estimation problem 
for the inter-jumping times. First, we define the estimator f”'{x,t) of f{x,t) by. 


V {x,t) e F, f^{x,t) 


F"{x,t) 

v^{x) 


with the usual convention 0/0 = 0. We have the following result of convergence. 


Corollary 3.8 For any couple {x,t) £ F such that Voc,{x)f{x,t) > 0, for any {a, [3) G A and (vo,wo) 
satisfying (14), we have 


f"{x,t) ^ f{x,t) 


and 


l-ad-& 


(/" {x, t) - fix, t)^ AT ^0, 


X? rj fjx, t) 

(1 -I- Qd -I- /3)j/oo(x) 


Proof. This result is a direct application of Theorem 3.3 and Slutsky’s lemma. □ 

Another feature of interest for the inter-jumping times is the survival function G. One may estimate this 
quantity by, 

V(x,f) G F, G"(*,t) = 

with the convention 0/0 = 0. Some properties of convergence are stated in the following corollary. 


Corollary 3.9 For any couple {x,t) G F such that Uoo{x)fix,t) > 0, for any a such that ad < 1 and 
a{d -f 2) > 1, and any vq, we have 


G"{x,t) ^ Gix,t) 


and 


(d"ix,t) - Gix,t)'^ -^Ar(o, 


rj Gjx, t) 

(1 -|- ad) 

^OO ix) 


Proof. The almost sure convergence is a direct application of Theorem 3.3. The central limit theorem 
is a consequence of the asymptotic normality established in Remark 3.6 and Slutsky’s lemma. □ 


3.2 Optimal estimation of the jump rate 

We propose to estimate the composed function A o $ by the ratio A o $ defined by, 

y{x,t)£F, Xo^ {x,t) = C^, (16) 

g"{x,t) 

again with the convention 0/0 = 0. Pointwise convergence and asymptotic normality are again a conse¬ 
quence of Theorem 3.3. 
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Corollary 3.10 For any couple {x,t) € F such that Vao{x)f{x,t) > 0, for any {a, fi) € A and (vo,wo) 
satisfying (14), we have 


A o $ {x, t) —4 A o 4>(a:, t), 


and 


l-ad-0 


^A o $ (a;, i) — A o <3>(a;, t)^ J\f ^0, 


T~i T-J A o $(a:, t) 

(1 + ad + P)vac,{x)G{x,t) 


Proof. This result is a direct application of Theorem 3.3 and Slutsky’s lemma together with (6), (8) 
and (9). □ 

In all the sequel, we focus on the estimation of the rate A (and not on the composed function A o $) 
from the estimate Ao$ defined in (16). In particular, for some fixed value x G E, we introduce a 
class of estimators ^•^(x) of A(x) indexed by the elements ^ of the curve Cx described by the reverse 
flow 4?(r, —t), t > 0, and we propose to choose an estimate from this class in an optimal way. Before 
proceeding further we define the notation Cx as 


Cx = {^(x,—t) : 0<t<t (a;)}. 

By definition (1) of t~{x) we have Cx C E. In addition, for any ^ € Cx, we define Tx{^) as the unique 
time satisfying 

= ^€Cx. 

Thus we have the following trivial result, 


V^£C., Ao4>(C,r40) = A(a;). (17) 

As a consequence, we propose to define a class of estimators of A(a;) by 

A"(a;) = {A^(a;) = ^ ^ Cx} , (18) 

where A o $ estimates A o <0 (see Corollary 3.10) and has already been defined in (16). By virtue of 

Corollary 3.10 and only using (17), one has, as n goes to infinity, 

V^GCx, X^{x) \{x) and n 2 ^ ^A^(a;) — A(a;)^ Af (O, CT|,,:,_^(a:)) , 

where the asymptotic variance is given by 

with Kxi^) = nooiOG{^,rx{^)). (19) 

In this paper, we choose to approximate A(x) by the element A^(x) G A”(x) minimizing the asymptotic 
variance a^^^ p{x). In other words, our optimal estimator of the jump rate A(x) is obtained as 

A*(x) = A^, (x), where = argmax Kx{(,). 

S.&Cx 


Remark 3.11 A good criterion should be to maximize the invariant measure ^ G Cx- Indeed, 

if ^oo(^) is large, a large frequency of post-jump locations around ^ may be available in the dataset. 
Nevertheless, roughly speaking, the quantities of interest /(C)''4(0) G{f,,Tx{^)) are well estimated if 

a large number of post-jump locations are around C together with inter-jumping times around Tx{^). This 
naive criterion may he corrected by including the quality of the estimation at time rx{^), that is to say, 
by maximizing the product Vaa{0G{^,Tx{^))■ We also refer the reader to the simulation study presented 
in Section 4 and more precisely to Figure 8. 

The criterion Kx{^), C Cai) is generally uncomputable from the known features of the PDMP and 
thus remains to be estimated. In light of Theorem 3.3 and by definition (19) of Hxi^), we naturally 
propose to approximate this quantity by, 

VCGC., K"(C)=e"(C,r4C)). (20) 
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As a consequence, we propose to estimate the jump rate A(a;) by its statistical approximation in A"(x) 
maximizing the estimated criterion More precisely, 

\f{x) = X~n{x), where = argmax (21) 

The high oscillations or alternatively the high smoothness of a kernel estimator with a ill-chosen band¬ 
width suggest that the choice of the parameter a appearing in k” (see eq. (12) and (20)) is crucial at 
this maximization step. 

3.3 How to choose bandwidth parameters a and (37 

This part is devoted to the choice of the bandwidth parameters a and p. The criteria that we introduce 
in this part to choose these features are defined as line integrals along the curve C^- As a consequence, 
we need to ensure that this kind of quantity is well defined in our setting. 

Assumption 3.12 For any starting condition ^ G E, the reverse flow t i—<E>(^, —t) defines a change of 
variable, that is to say, is a diffeomorphic mapping. 

It is common in the literature to minimize the Integrated Square Error (ISE) to choose the optimal 
bandwidths of a kernel estimator, in particular from dependent data: one may refer the reader to 
[23, 26, 27]. Another classical solution is to investigate the behavior of the Mean Integrated Square Error 
(MISE). In the framework of Gaussian dependent data, the authors of [15] have shown that the optimal 
bandwidths obtained by minimizing the ISE and the MISE are very close if the dependence is of short 
range, that is to say, if the covariance function is integrable (see [15, Theorem 2.2]). It should be noted 
that a geometric ergodic Markov chain satisfies this kind of condition (see [29, Theorem 16.1.5]). 

Let us recall that a appears at first in the computation of the estimated criterion k" which we need 
to maximize along the curve C^. Indeed k" is computed (20) from the estimate Q'^ which implicitly 
depends on a. We propose to choose the bandwidth parameter a by minimizing the ISE associated with 
Kj, and defined by 

ISE::(a) = ^ («"(?)-K.(C))'d^ 

= J Kxiifdi + e«(a), (22) 

where the function is given by 

f g"(C,r40)K40de (23) 

Jc,, 

One may remark that here the ISE is unusually computed along a curve of interest. In (22) the dependency 
on a only holds through the function e”. As a consequence the optimal parameter a minimizing this 
stochastic function (23) also minimizes the ISE. 

The function e" is generally not computable since the unknown quantity appears in its definition. 
As a consequence, we propose to estimate e" by cross-validation which is a popular technique for se¬ 
lecting the bandwidth that minimizes the ISE. The authors would like to highlight that cross-validation 
involves here two main difficulties. First the estimators are computed from dependent data which are not 
identically distributed. In addition, there is almost surely no data on the set of integration Ca, whenever 
the dimension d is larger than 2. That is why we propose a specihc procedure adapted to this framework. 

Before dehning our cross-validation estimate of e™, we need to introduce the quantities ^x,p and 
First IHa; denotes the hyperplane orthogonal to Cx at x, that is, 

= : y - a; T Vt3>(2;, 0)} . 

In addition, for any p > 0, we introduce the notation Dx,p for 

IDa:,p = Bd{x,p) n IHa,. 
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Furthermore ^x,p denotes the tube around Cx with radius p, 


~^X,P — U 


( 24 ) 


Finally, for any ^ dx{C) denotes the unique time such that ^{^,9x{^)) £ '^x,p- In particular, 

V^GC,„, = (25) 

It should be noted that, for p small enough, lx,p C E. 

We focus now on a cross-validation method for estimating the quantity e". The estimate eof e" 
is defined from the observation of the embedded Markov chain (Z*,, Sk+i)k>o of another PDMP {Xt)t>o 
(independent on the first one (Xt)t>o and distributed according to the same parameters), by 


^"’'>(a)= / - ^^Wz,,04Z,)) lT.,p(^0l(e.(z,),+oc)(5.+i), 

n-TT 2 nd-^ ^ ^ 


(26) 


where F denotes as usually the Euler function. Some regularity conditions are necessary to investigate 
the asymptotic behavior of the cross-validation estimate in Proposition 3.14. 

Assumptions 3.13 

• The sup-norm ||:^oo||oo is finite and v^o is Lipschitz, 

y X, y £ E, \voa(,x) - i'oc,{y)\ < [i'oa\Lip\x - y\. 

• The deterministic exit time for the reverse flow t~ is Lipschitz, 

yx, y £ E, \t~{x) -t~{y)\ < [t~]Lip\x - y\. 

• The flow 'F is Lipschitz, 

yx,y £ E, yt£R, |4>(a;,t) - ^{y,t)\ < [<E>]Lip|a: - y\. 

• The sup-norm llVt^Hoo is finite and Vt*!* is Lipschitz, 

yX, y £ E, Vt G R, |Vt$(a;,t) - Vt$(j/,t)| < \Vt^]Lip\x - y\. 


Proposition 3.14 Conditionally to a(Zo, Si, ... ,Z„_i, Sn), we have 

Jim = £"( 0 ) a.s. 

n —^00 

p—>-0 

Proof. The proof is stated in Appendix C. □ 

By virtue of Proposition 3.14, one may obtain an estimate of the optimal bandwidth parameter a 
arising in Hf, by minimizing the quantity ej’”’'’ for some small enough p and large enough n. In addition, 
by (16), (18) and (20), the quantity kJ also appears in the calculus of the estimator A^(x). In particular, 
the same choice of a may be done for computing the denominator k.x{0 oI -^^(s;)- As a consequence it 
remains to choose in an optimal way the bandwidth parameters a and fi arising in the formula (16) of 
the numerator E'^{fi,rx(fi)') of Aj(a:) (18). 

In a similar way as before, we propose to choose ol and /3 by minimizing the ISE associated with 
.^"(■i'ra;(-)) and computed along the curve Cx, 

ISE^(a,J) = ^ (.^"(e,T4C))--^(C.r.(e)))"d^ 


= L 


Tx{i)) di + £]r(a, 15), 


(27) 
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where T" implicitly depends on a and /3, stands for Vrx{i)f and is given by 

e-Aa,0)= f ^"(C,r40)"dC-2 / ^"(C, ^4e))dC. 

As in the previous part we propose to estimate by cross-validation from the observation of the 
embedded chain {Zk, Sk+i)k>o of another PDMP {Xt)t>o- We define our estimate qJ- 

^ -F"(e,r40fde 

2p C rf-i I ^ ^ „ _ 

“ Z- Li d_i )+fa)(5fc-ri)- (28) 

np2TV 2 k=o 

The convergence of is investigated in Proposition 3.15. 

Proposition 3.15 Conditionally to a-(Zo, Si, ... ,Z„-i, Sn), we have 

lim P) = £^{a,P) a.s. 


n—^oo 
Pi,2^0 


Proof. The proof is similar to the demonstration of Proposition 3.14 stated in Appendix C. 


4 Simulation study 

In this section, we provide a self-contained presentation of the estimation procedure, as well as three 
application scenarios on both simulated and real datasets. 


4.1 Estimation algorithm 

The sequel is devoted to the self-contained presentation of the estimation procedure provided in this 
paper. Precisely, we are interested in the estimation of X{x) for some x £ E <Z R'*. 

4.1.1 Preliminary computations 

These preliminary computations only require to manipulate the flow <1? and the state space E. 

• Compute the curve Cx = —t) : t > 0} n i5. 

• Choose p > 0 and compute 

^x,p = Bd{x,p)c{y : y - x ± Vt${x,0)}, 

~^X,p = [J Cy. 

y^^x,p 

• For any ^ € Tx,p compute as the unique solution of 

HMO) e D.,p. 

The mapping Tj, is only defined on Cx C ~^x,p by Tx{^) = 9x[p). 


4.1.2 Preliminary estimates 

These preliminary computations require to choose kernel functions Ki and Kd on R and R'*, as well as 
two positive values vo and wo¬ 
rn For any couple {a,P), compute 




1 

nv^wo 


n—1 

i=0 



Ki 



• Compute 




1 

nug 


n—1 

i=0 



l{Si+i>t}- 
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4.1.3 Choice of bandwidth parameters by cross-validation 

From the observation of two independent embedded chains {Zi, Si+i)o<i<n-i and {Zi, Si+i)o<i<n-i, 
with n n, one determines the optimal bandwidth parameters appearing in the preceding estimates. In 
practice, from only one trajectory of the underlying PDMP, one may divide the data into two categories: 
the largest one is used for the estimation, while the cross-validation step relies on the other one. 

• Compute 


: argmax 

Q :>0 


L 





,iZk) 1 


(Sa:{Zfc), + Oo) 


(Sk+i). 


• Choose p 2 > 0 and compute 


L 


{a ) = argmax / .F'a,/ 3 (C'rx(C)) d^ 

c«,/3>0 


2 p C _|_ ij ^ _ 

3 d-1 J^a,pyZk,dx{Zk)j + 


np2n 2 


4.1.4 Estimation of A(a:) 

Finally, we compute the best estimate of A(a;) from all the preceding computations. 

• Compute 

C = argmaxg"s(CT:,(0)- 

• From {a^,l3^), cfi and compute 

Remark 4.1 The time complexity of this algorithm depends on several parameters, namely the number 
n of observed jumps, the number h of observed data for the cross-validation steps, and also the numbers 
N^, Na and Ng for discretizing the state spaces of C o ft at each maximization procedure. It 
is easy to see that the cross-validation step is the most complex but remains polynomial, precisely in 
OinNlNpiN.+n)). 


4.2 TCP-like process 

The application on which we focus in this part is a variant of the famous TCP window size process 
appearing in the modeling of the Transmission Control Protocol used for data transmission over the 
Internet and presented in [12]. This protocol has been designed to adapt to the traffic conditions of the 
network: for a connection, the maximum number of packets that can be sent is given by a random variable 
called the congestion window size. At each time step, if all the packets are successfully transmitted, then 
one tries to transmit one more packet until a congestion appears. 

The model presented in this part is two dimensional. For the sake of clarity, we will use the following 
notation: for any x G R^, Xi and X 2 denote the components of x. We consider a PDMP {Xt)t>o evolving 
on the state space E = (0,1)^. The deterministic part of the model is defined from the flow given by, 

V® £ R^, Vt £ R, 4>(a;, t) = (xi -I- t, X 2 ). 

The jump rate A is defined by, 

Va: £ R^, X{x) = x\ X 2 . 

The transition kernel Q is defined by, 

'ix G E, y Ag B{E), Q{x,A)<x [ u{l - v{l - v) dudv. 

J A 
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Starting from {x-i,X 2 ), the process evolves in the unit square, always to the right, until a jnmp appears 
either when the motion hits the boundary {(^ 1 , ^ 2 ) : Ci — 1} or with the non homogeneous rate a;i+ 2 : 2 +t 
before, that is according to a Weibull distribution. The two components of the post-jump location are 
independent and both governed by a Beta-distribution, in such a way that the process tends with a high 
probability to jump to the left of the location just before the jump. One obtains a TCP-like process (see 
Figure 2) for which the second dimension models the quality of the network (upper the second component 
is, higher the probability of a congestion is). 

[ Figure 2 approximately here ] 

The asymptotic behavior of the process may be represented by the invariant distribution Uao of the 
post-jump locations. Since this quantity is unknown, we propose to show in Figure 3 its estimate 
defined by (13) and computed from n = 20 000 observed jumps. 

[ Figure 3 approximately here ] 

We present here all the procedure for estimating the jump rate A at the location x = (0.75,0.5) for 
which the quality of the network is average. In this context, the class of estimators of A(a:) is indexed by 
the elements ^ = (0,0.75] x {0.5}. It should be noted that the invariant distribution at x is quite 

low (see Figure 3). Nevertheless our method is expected to work pretty well even in this unfavorable 
framework. 

In this simulation study, we assume that we observe the embedded Markov chain [Zn, Sn+i)n>o until 
the 10 000*** jump. The cross-validation procedure is computed from an additional chain, independent 
on the hrst one, and observed until the 1 000*** jump. When boxplots are presented, they have been 
computed over 100 replicates. 

We begin with the choice of the bandwidth parameter a appearing in k"(^) = The 

cross-validation procedure relies on the minimization of the estimate £'"’"’^(a) which depends on the 
positive parameter p. We present in Figure 4 this quantity as a function of a and from different values 
of p. Fortunately, this new parameter seems to have little influence over the behavior of the estimation 
of the ISE along Cx- 

[ Figure 4 approximately here ] 

Now, from this a (denoted in the sequel a®), we maximize the estimated criterion along the 

curve Cx (see Figure 5): we obtain the optimal point Q at which we will compute our estimator of A( 2 ;). 
The crucial role of a at this maximization step is illustrated in Figure 6. 

[ Figure 6 approximately here ] 

We continue with the choice of the couple {a, j3) implicitly appearing in the estimator .F’*(^, rj,(^)). 
The optimal parameters (denoted in the sequel {oL^,j3^)) are obtained by minimizing the estimate 
of the related ISE (see Figure 7). 

[ Figure 7 approximately here ] 

We compute the estimators A^ (a;) for different values oi ^ £ Cx and with the optimal bandwidths q® 
and The related boxplots are presented in Figure 8. The procedure makes us able to choose 

the best index Q which most of the time corresponds with the estimate with least bias and variance. 
This proves the strong interest of the estimation algorithm developed in this paper. 

[ Figure 8 approximately here ] 


4.3 Bacterial motility 

We present here a model of bacterial motility. Most motile bacteria move by the use of a flagellum or 
several flagella. The bacteria moves in the direction of the flagellum. A flagellum behaves like a rotary 
motor. Periodically, the flagellum changes its direction and results in reorientation of the bacteria. This 
allows bacteria to change direction. Bacteria can sense nutrients and move towards them. Additionally, 
they can move in response to temperature, light, etc. If the bacteria is in a favorable environment, 
the frequency of changes in direction is low. This intelligent behavior is allowed by the fact that the 
jump rate of the direction depends on the environment. For example, [30, 31, 34] propose models for 
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the trajectory of the bacteria E. Coli. In particular, the author of [30] uses PDMP’s to describe the 
movement of bacteria under the influence of an external attractive chemical signal. 

We present here a variant of the model presented in [30]. The path of the bacteria is described by 
a PDMP {Xt)t>o evolving in the three-dimensional space state E — D x [0, 27r) with D the unit disk 
of R^. For the sake of clarity, we will use the following notation: for any x ^ , xi, X 2 and 0:3 denote 

the components of x. In our case, (xi,X 2 ) is the position of the bacteria in the unit disk, and X 3 is the 
direction of its flagellum. 

The bacteria moves in the direction of its flagellum. In other words, the flow is given by, 

Va: £ R^, Vt £ R, ^{x,t) = {xi + tcosxa, X 2 -I- t sin* 3 , * 3 ). 

The bacteria changes its direction according to a jump rate function X{x) = X{x\,X 2 ) which depends 
on the environment only through the position (a;i, 2 : 2 ) and not on the current angle. The jump rate 
describes the interaction between the bacteria and its environment. 

When the bacteria changes its direction, the new direction is chosen preferentially in the direction 
of a more favorable environment. The transition kernel Q{x, •) models this change of direction. In our 
case, we suppose the bacteria has no information a priori on the quality of the environment around itself. 
Thus, Q can be defined by, 

y X € E, MA £ I?([0, 27r[), Q(x, {xi,X 2 ) x -4) = f dw. 

Ja 

Starting from the position {xi,X 2 ) in D, the bacteria evolves in the direction X 3 £ [0,27r), until a jump 
appears either when the bacteria hits the boundary dD of its environment or with the rate A (‘I>( 2 ;, t)) = 
A(a:i -I-1cos* 3 , 0:2 -|-tsin 2 ; 3 ) before. The post-jump direction is next randomly chosen in [0,27r), and the 
bacteria continues its path in this new direction, and so on. A possible path of this model is presented 
in Figure 9. 

[ Figure 9 approximately here ] 

In this simulation study, the jump rate is taken constant equal to A = 1. In other words, there is 
no interaction between the bacteria and its environment. Of course, this assumption is not taken into 
account in the estimation. Our goal is to estimate the jump rate A in different points of E in order to 
check if it depends on the position. Actually, this provides an estimate of the influence of a likely external 
attractive signal on the bacteria. 

The estimation algorithm presented in this paper yields an estimate of the jump rate at the three- 
dimensional state {xi,X2,0). Nevertheless, one should take into account that the jump rate do not 
depend on the angle 9. We explain here the procedure to reduce the dimension and estimate the jump 
rate at the location (a;i, 2 : 2 ) £ D. 

• Let 9 £ [0, 27r) be a fixed angle. Estimate the best index ^t{9) and the jump rate X~„^g^(xB) at the 
location xg = (xi,X2,9). 

• For any 9, X~^^g^(xg) estimates A(a:i,a: 2 ). We reproduce the preceding step for any 9 £ [0,27r). We 
obtain both a connected trajectory around (a;i,X 2 ) 

£.„., = {C(9):9elO,27r)}, (29) 

and a class of optimal estimators 

: 0£ [0,2^)} = {A^ : ^ 

indexed by elements of £xi, 2 : 2 - 

• As the jump rate does not depend on the direction 9, we aggregate the estimators in to 

obtain the following jump rate estimate 

1 

X{xi,X 2 ) = —J A|„(g)(x 0 )d 6 l. (30) 
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We investigate the estimation of the jump rate at the following 9 target points, 

(* 1 , X 2 ) G {(0,0), (-0.5,0), (-0.5,0.5), (-0.5, -0.5), (0,0.5), (0, -0.5), (0.5, 0), (0.5, 0.5), (0.5, -0.5)}, 

from a trajectory of n = 100 000 changes of direction. In the sequel, the aggregated estimator (30) is 
approximated from a discretization of the interval [0, 2tt) with step tt/S. For each of the 9 target points, 
we present in Figure 10 the boxplot of the estimates for 0 in the discretization grid. 

[ Figure 10 approximately here ] 

As shown in Figure 11, the trajectories £xi,x2 defined in (29) are very similar and close to the 
boundary dD. We explain this fact as follows. For each target state x = {xi,X2,0), the optimal point 
(6») maximizes k"{^) = along 

Cx = {{xi — tcos 6 ,X 2 — ts'm9,6) : t>0}nD, 

where Ta;(^) satisfies xicos0 = and * 2 -'r 2 :( 5 ) sinS = ^ 2 - Recall that ra;(^)) is an estimator 

of lixiO — ^oo(C)G'(^, Xx{^)). On the first hand, G(^, t) = exp(—t) for t small enough, then the only point 
that maximizes G(^, T"a;(^)) is ^ = a; itself. On the other hand, there is an important number of data that 
are very close to the boundary dD (due to the fact that the bacteria jumps when it hits the boundary 
dD). Then the maximization of G”'{^,Tx{^)) is a compromise between ^ = x and ^ near the boundary. 

Finally, the estimates A(a:i,X 2 ) at the 9 target points are close together, and close to the true value 
A = 1 with a small variance, which is a good indicator of a constant jump rate (see Figure 11). 

[ Figure 11 approximately here ] 

Now, we investigate the effect of the size n of the dataset on the quality of the estimation. We 
compare results of simulations from n — 100 000, n = 50 000 and n = 20 000 observed data. When 
n = 20 000, the trajectories €xi,x2 are not all the time close to the boundary. It is due to the fact that, in 
this case, some areas of the state space are not well-covered by the observed post-jump locations. Then, 
if the target point is x = {xi,X 2 , 9 ), the optimal point ^"{6) € Cx is obtained at the location x itself (see 
Figure 12). In addition, the estimated curves £xi,x2 obtained from n = 50 000 data are very close to the 
ones computed from n = 100 000 observed jumps. 

[ Figure 12 approximately here ] 

Of course, the larger the number of data n is, the better the estimations (see Figures 13 and 14). In 
particular, it is difficult to conclude that the jump rate is constant from only n = 20 000 observed jumps 
because the variance of the estimates is too large. Nevertheless, a dataset of size 50 000 may be sufficient 
to conclude that the jump rate does not depend on the location. 

[ Figures 13 and 14 approximately here ] 


4.4 Fatigue crack propagation 


Fatigue crack propagation is a stochastic phenomenon due to the inherent uncertainties originating 
from material properties and cyclic mechanical loads. As a consequence, stochastic processes offer an 
appropriate framework for modeling crack propagation. Fatigue life may be divided into crack initiation 
and two crack growth periods, namely the linear or stable regime described by Paris’ equation 

|^=G(AA(a))’", 

and the acceleration or unstable regime modeled by Forman’s equation 

da _ C{AK{a))^ 

diV ^ (1 - R)Kc - AK{a)’ 

where a is the crack length, N is the time measured in number of loading cycles, R is the stress ratio, 
Kc represents the maximal value of the stress intensify factor required to induce failure and AK{a) is 
the range of the stress intensity factor given by 


AK{a) — Act cos 
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where to is the size of the test specimen and Act is the stress range. In addition, m and C are two 
unknown material parameters. In this context, regime-switching models have been proposed to analyze 
fatigue crack growth data [7, 14]. The authors of [7] propose to estimate the transition time between 
Paris’ and Forman’s regimes assuming that the crack propagation follows the trajectory of some piecewise- 
deterministic Markov process. The estimation procedure is performed on Virkler’s dataset [35] in which 62 
identical centre-cracked aluminium alloy specimens (oj = 152 mm) were tested under constant amplitude 
loading Act = 48.28 MPa at a stress ratio R = 0.2. The number of loading cycles for the crack tip to 
advance a predetermined increment Aa was recorded from an initial crack length of no = 9 mm to a final 
length of 49.8 mm. 68 crack growth histories were obtained from these tests (see Figure 15). 

[ Figure 15 approximately here ] 

In this paper, we assume that the crack growth propagation follows Paris’ equation with random 
parameters m and C before switching to Forman’s equation with another set of parameters m and C 
distributed according to some transition measure. The transition occurs at a random time Ti given by 
its survival function 


P(ri > t \ ao = 9 mm, m, C) = exp 




where is fh® deterministic flow of Paris’ equation with parameters m and C. In the sequel, 

will be computed from the Runge-Kutta method because there is no explicit solution to this differential 
equation. We propose to estimate the jump rate of this stochastic model from real crack growth data. 
It should be noted that the jump times are not directly observed in Virkler’s dataset but have been 
estimated in [7]. The jump rate is of significant importance in understanding the transition between the 
stable and unstable regimes of crack propagation, and has never been estimated in the framework of 
PDMP’s. We emphasize that Virkler’s dataset is composed of 68 independent experiments and thus does 
not directly follow the theoretical framework developed in this manuscript. Nevertheless, our estimation 
procedure is expected to perform pretty well in this more favorable context (independent curves instead 
of only one Markov path). 

Material parameters m and C are unobserved but have also been estimated from Virkler’s experiments 
in [7]. It is well-known in the literature that there is a strong linear relationship between these features. 
This characteristic has again been highlighted in [7] and is shown in Figure 16. In order to reduce 
the dimension and thus simplify the model, we propose to parametrize Paris’ equation with only one 
parameter, say m, while log(C) is obtained from this linear relationship. 

] Figure 16 approximately here ] 

We focus on the estimation of A(a) for some crack length a. For each parameter m of Paris’ equation 
there exists a unique deterministic time Tm{a) such that Paris’ flow reaches a at time Tm(a) starting 
from ao = 9 mm. Criterion Ka(m) and its maximum have been estimated for different values of a (see 
Figure 17). In particular, one may observe that the larger the target length a is, the larger the optimal 
parameter m. 

[ Figure 17 approximately here ] 

For each crack length a, we estimate the jump rate X{a) from the parameter m maximizing the 
estimated criterion Ka{m). We obtain the estimated function displayed in Figure 18. This curve is 
increasing as expected and describes the transition rate from the stable region of propagation towards 
the acceleration regime that leads to fracture. This makes us able to detect conditions of crack growth 
instability and could be used to predict the critical length in fatigue crack propagation with a given level 
of confidence. 


] Figure 18 approximately here ] 


A Ergodicity and invariant measures 

In this section we present some preliminary and technical results about the invariant distributions of the 
underlying Markov chains. We begin with some properties of the Markov chain {Zn)n>Q of the post-jump 
locations of the PDMP of interest. In particular, under Assumption 2.1, one may state the following 
result. 
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Proposition A.l We have the following statements: 

• {Z„)n>o is Vao-irreducible, positive Harris-recurrent and aperiodic. 

• i/oo is the unique invariant distribution of {Zn)n>o. 

Proof. The proof is similar to the demonstration of Proposition 4.2 of [4]. □ 

These properties make us able to apply the law of large numbers to the Markov chain {Zn.)n.>o (see [29, 
Theorem 17.1.7]). Now we propose to focus on the sequence {Zn, Sn+i)n>o which also forms a Markov 
process whose transition kernel is given by TZ in (11). Let us recall that pn denotes the distribution on 
the state space F defined by (10) of the couple {Zn, Sn+i), n > 0. We define the measure poo by, 

y A G B{R‘^), y t > 0, pao{A X {t,-\-oo)) = I Voo{dx) S{x,{t,-\-oo)), (31) 

J A 

where the conditional distribution S is given by (6). In particular, pao is the unique invariant distribution 
of {Zn, Sn+l'ln^Q. 

Proposition A.2 We have the following statements: 

• For any initial distribution po = {x,t) G F, we have 

lim ||/i„ —/ioojjTV = 0. 


. {Zn,Sn +i)ti>o I’S Poo 'irreducible, positive Harris recurrent and aperiodic. 

• Poo IS the unique invariant distribution of {Zn, Sn+i)n>o- 
Proof. The proof is similar to the demonstrations of Lemma 4.4 and Proposition 4.5 of [4]. □ 

We are able again to apply the law of large numbers to the Markov chain {Zn, Sn+i)n>o by virtue of this 
result. It should be noted that the measures pn, n > 0, and poo do not share with Vn, n > 0, and noo 
the property of absolute continuity with respect to the Lebesgue measure presented in Remark 2.2. 

Remark A. 3 For any n, pn and thus the limit poo admit a density only on the interior F of the state 
space F because of the expression (7) of S and thanks to Remark 2.2. 

Finally, we would like to highlight that the link between the measures Voo and may be expressed in 
another way than (31). Indeed, for any x G E, we have 

/ Q{^{y,s),x)poo{dy X ds) = noo{x), (32) 

J F 

by the expression (5) of the transition kernel V of {Zn)n>o and because i/ooP = i^tx 3 by virtue of Propo¬ 
sition A.l. The formula (32) of i/oo will be useful in our investigations. 


B Proof of Theorem 3.3 

In all this section we use the notation for the cr-algebra a{Zo, Si, ... , Zn, Sn+i), n > 0. In addition, 
the classical symbols ~, o and O must be understood to hold here in the almost sure sense. 


B.l Sketch of the proof 

The proof of Theorem 3.3 relies on the following decomposition. For any integer n > 1, we have 


■ F"{x,t) ' 


Voo{x)f{x,t) 

§"{x,t) 

- 

noo{x)G{x,t) 

. i^S>(a:) . 


noo{x) 


(33) 


where the sequence {Mn)n>o is a (F„)n>o-vector martingale defined by (37) and studied in Appendix 
B.3, and TZn is a remainder term defined by (34) and studied in Appendix B.2. 
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In Appendix B.2, we establish that the remainder term almost surely goes to 0 when n tends to 
infinity. This is a first step to show the almost sure convergence presented in Theorem 3.3. In addition 
we investigate the rate of convergence of the remainder term in Appendix B.2 under an additional 
Lipschitz mixing condition stated in Assumption 2.4. This is enough to prove the asymptotic normality 
given in Theorem 3.3. 

The rest of the proof deals with the study of the martingale term. We prove in Appendix B.3 that the 
process (A4„)„>o is a vector martingale. We investigate its asymptotic behavior by studying its square 
variation process ((M}n)n>o in Appendix B.3. Thanks to these results we state the law of large numbers 
and the central limit theorem for {Mn)n>o in Appendix B.3. 

Finally the almost sure convergence presented in Theorem 3.3 is a direct application of (33) and (35) 
together with (51), while the asymptotic normality is obtained from (33), (36) and (54). 


B.2 Remainder term 

This part of the paper is devoted to the asymptotic study of the remainder term sequence {TZn)n>i 
appearing in (33). 

B.2.1 Definition of the remainder term 

For any n, one may write TZn = , where each of the components is defined by 

1 Jl<fc<3 


= E ’ = E “d = E 


(34) 


The five terms that define are given by 


^( 1 , 0 ) 

^( 1 , 1 ) 


nvgwo 

n 


Zo — x\ f Si — t 
1 


Vo 


Wo 


R\ 


( 1 , 2 ) _ 


R\ 


(1,3) _ 


E [ + vwj)Kd{u)Ki{v) {Q (<1> {Zj-i, Sj) ,x + uvj)-Q{^ ,x)) dudv, 

Jr‘‘xr 

n — 1 p 

E Q ('!’ {.Zj-i,Sj) , x) / [f{x + uvj,t + vwj) - f{x,t + vwj)] IKd(u)IKi(u) dudu, 

ilR^xR 

E Q ('!’(^i-ii S'i), ai) / [f{x,t + vwj)-f{x,t)]Kd{u)Ki{v)dudv, 

Jr'^xR 


R\ 


(1,4) _ 


i=i 

n — 1 


1 ^ 

-Y,Q{^{Z,.i,S,),x)-v^{x) 


j = ^ 


fix,t). 


( 2 ) 

The four terms that define R„ ' are given by 

Ir 


r(2,0) 

^( 2 , 1 ) 

l-Tl. 




n vl 

n — 1 


Vo 




r: 


( 2 , 2 ) _ 


1 “_^ r 

- E / {Zj-i,Sj) ,x + uvj) - Q ('F {Zj-i,Sj) ,x)] G{x + uVj,t)Kd{u) du, 

1 r 

-E Qi^ / lG{x + uvj,t) - G{x,t)]Kd{u) du, 

n J^d 


R. 


(2,3) _ 


f = l 

n — 1 


-Y.Q{^{Zi.i,Si),x)-v^{x) 

T) ‘ ^ 


j=i 


Gix, t). 
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f 3^ 

Finally, the three terms that define Rn ' are given by 

r(3.0) 


1 


R 


R 


Vo 


(3.1) _ 


(3.2) _ 


1 r 

1 


1 = 1 


B.2.2 Almost sure convergence 

We only investigate the first component (i?£^E>i since the other terms may be treated with similar 
arguments. First, it is obvious that the sequence almost surely tends to 0. In addition, 

the term converges to 0 by virtue of the ergodic theorem applied to the Markov chain 

(Z„, Sn+i)n>o (thanks to Proposition A. 2 ) together with (32). Because both the functions Q and / are 
Lipschitz and bounded (see Assumptions 2.3), we have 

I [,(1,1)1 ^ [21i*pll/lioo f I in/ [ 

.Rn < > Vj U IKd(li)du, 

II ^ 

|^(1.2)| < IIQIUl/l^ip g r 

II « ^ iR3 

/ |w|IKi(u)du, 

II ^ 

which are some Cesaro means because both the sequences (v„)„>o and (w„)„>o tend to 0. Thus, 

7^n ^ O 3 . (35) 


B.2.3 Rate of convergence 


Under the first item of Assumptions 2.4, the Markov chain (Zn)n>o and thus the two-dimensional process 
(Zn, S'n+i)n>o Satisfy the contraction assumption given in Theorem 6.3.17 of [21]. By applying this 
theorem to the function <fi = Q($(-, ■), ■) € Li(ri,r 2 ) (see second item of Assumptions 2.4), we obtain 
that, for any 7 G (0,1) and x £ E, 

72 — 1 

n-^ E {Zj-i,Sj),x) - u^{x)] ^ 0, 

1=1 


because 2(ri -I- 72 ) < cii. Therefore, we have for any couple {k,j) in the set {(1,4), (2, 3), (3, 2)}, 


l-ad-0 

n 2 


R 


(fc.!) 


0 . 


The same result is obvious for the other couples {k,j) with the condition ad -\- fd + 2min(a,/3) > 1. 
Finally we obtain, when n goes to infinity, 

n 5 7^„ ^ O 3 . (36) 


B.3 Martingale term 

B.3.1 (Afn)n>o is a vector martingale 

Let (Mn)n>o be the three-dimensional process defined for any n by j where each 

^ L Jl<fe<3 

component is defined as 

n 

mW = EE’ ^ {1,2,3}, (37) 

1=1 
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and the terms and are given by 


A' 


( 1 ) _ 




Sj+i — t 


B 


( 1 ) 


,{ 2 ) 


JIR'^XR 


Q{^{Zj-i,Sj),x + uvj)f{x + uvj,t + vwj) Kdiu)Ki{v) \d{du) dn, 


Zj X 




B 


( 2 ) _ 


/ Q{^{Zj-\,Sj), X + uvj)G{x + uvj,t) IKd(M) Ad(du), 

jRd 


A' 


(3) 




B 


(3) _ 


[ Q{^{Zj-i,Sj), X + uvj) Kd(?i) Ad(du). 


We claim that the process {Mn.)n.>o is a (F„)„>o-martingale. The keystone is to show that, for any n 
and fc, we have 

E [aI") I F„_i] = 

The proof presents no particular difficulty except for the first component fc = 1 for which we provide 
some details. Let us recall that TZ denotes the transition kernel of the Markov chain [Zn, Sn+i)n>o (see 
Subsection 2.1). We have 


E I F„_i] 


iR'^xRj- \ / \ Wn J 


V-nWn jRtixR^. 

1 


V-nWn jRtixR 


Vn 

f 

JRdxR, \ V„ J \ Wn J 


El 

El 


Wr 

S — t' 


Tl{{Zn-\,Sn)A^ X ds) 
Q(<f>(Z„_i,5„),d05(C,ds), (38) 


with (11). Let us recall that the bandwidth sequence (v„)n>o is decreasing. Together with the third 
item of Assumptions 3.1, we obtain 


supp Ed 


With similar arguments, we obtain 


C supp Ed 


Vo 


C Bd(x,voS). 


supp El 


• - t 
Wn 


C supp El 


■ - 1 
Wo 


C Bi{t,woS) C — 00 , inf 

' ieBaix,voS) 


(39) 


with the condition (14) on the couple (uo,uio). By (38), (39) together with the expression (7) of 5 we 
obtain 


L I J Vj^Wn J^a 


Ed 


^-x 


IRdxR^ 

One may conclude by a change of variable and (3). 


El 


s — t 


Q(4>(^„_i,5„),dO/(C,s)ds. 


B.3.2 Predictable square variation process 

The asymptotic behavior of the martingale {Mn)n>o may be investigated by studying its predictable 
square variation process that we denote as usual {{M)n)n>o- At any time n, {M}n is the 3x3 symmetric 
matrix defined by 


{M}n = ^ [e I F,_i] - 


j<3 


Now, we calculate each coefficient of this matrix by beginning with the diagonal terms and the first of 
these (A4}n’^\ One has 
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where the terms and are defined by 

[ Qi^{Zj-i,Sj),x + uvj)f{x + uvj,t + vwj)Kl{u)Kl{v)dudv, 

f / Q($(^j-i,S'j),a; + MWj)/(2; + Miij,t+ nwj)IKd(u)IKi(?;)dudri^ . 

ViRt^xR / 

Since the functions Q and / are bounded (see Assumptions 2.3), we easily obtain that 


( 40 ) 


Let us introduce the additional notation, 

~ A ^ f 

T^'’ = - / Q{^{Zj-i,Sj),x + uvj)f{x + uvj,t + vwj)¥.l{u)Kl{v)dudv, 

^ ^ iR<*xR 

1 ^ f 

= - VQ($(Z,_i,5,),a:)/(2;,t) / Kl(u)K?(u)dudu 

Jr^xr 

2 2'^ 


by Remark 3.2 on the kernel functions Ki and Kd. Since both Q and / are Lipschitz and bounded (see 
Assumptions 2.3), together with suppKp C Bp{0p,5) (see Assumptions 3.1), the stochastic sequences 
(Ti^’ )n>0 and (ri^) )ti^o have the same limit thanks to 


n ^71 


^ m\\oo\fUp + [Q]L.AfU)Sr!r? ^ ^ ^ 

1=1 


(41) 


which is a Cesaro mean. In addition, by the almost sure ergodic theorem (see Proposition A.2) together 
with (32), we obtain 

f^'> ^ TlT'^Vaa{x)f{x,t). (42) 

By (41) and (42) and by virtue of Lemma 7.1.5 of [21], we obtain 

a.s. T]^TiU^{x)f{x,t) 

n°:d+f} i^ad + p 

Together with (40), as n goes to infinity, 

(-^)n a.«; 'rirfi^^{x)f{x,t) ^ 

^i+od+/3 i + ad + P ^ ^ 

The second diagonal term of the predictable variation process may be studied in a similar way as before. 
We obtain 


{M) 


( 2 . 2 ) 

n 


E 


1=1 


[ Q{^{Zj-i, Sj), X + uvj)G{x + uvj,t)¥,^^{u)d'. 

/ Q(4>(Zj_i, Sj),x + uvj)G{x + uvj,t)Kd{u)du 

jRd 


and, as n goes to infinity. 


{M)n'^'’ a.s. Tlvoa{x)G{x,t) _ 
p^l + txd ^ 1 + ad ^ 


(44) 
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The third and last diagonal term may also be investigated in the same way. We have 

\ Sj),x + uVj)Kaiu)du 

'j Jr‘‘ 


(3,3) _ 


E 

i=i 


and, when n goes to infinity, 


~ (^J + uvj)Kd{u)du 

a.s. rjv^ix) _ ^ 

^1 + ad ^ 1 , ^ j 


1 + ad 

Now we focns on the non diagonal terms. For any integer n, we have 


( 45 ) 


{M) 


( 1 , 2 ) _ 


E 


/IR^^xIR 


{M) 


(1,3) 


E 


Q{^{Zj-i, Sj),x + uvj)^f(x + uvj,t + vwj)G{x + uvj,t)Kd{u)Ki{v)du dw 

— ( / Q{^{Zj-i, Sj),x + uvj)f{x + uvj,t + vwj)Kd{u)Ki{v)dudv] 

vilR'ixIR / 

X Q(^{Zj-i, Sj),x + uvj)G{x + uvj,t)Kd{u)dii^ , (46) 

/ Q(4>(.Zj-i, Sj), a; + uvj)^f{x + uvj,t + u'u;j)IK^(u)IKi(ii)dwdu 

aiR^xIR 

-( / Q{${Zj-i, Sj),x + uvj)f{x + uvj,t + vwj)Kd{u)Ki{v)dudv^ 

\Jr‘‘xr / 

y^Q(4>(Z,_i,5,),a; + nn,)B<d(M)dM^ , (47) 


71 ' p 

(A4)i^’®' ^ E / Q(^(^i-i>‘S'j),a: + MVj)^G(® + un3,t)IKd(w)d« 

Q(4>(Zj-i, 5j), a; + uvj)G{x + uvj,t)Kd{u)dii^ 


X [ / Q{^{Zj-i,Sj),x + uvj)Kd{u)duj 


From (46), (47) and (48), and because Q and / are bounded (see Assumptions 2.3) together with the 
fact that the integral J IK^dAp is finite (see Remark 3.2), we easily obtain that, for any i ^ j, 

= 0(n). (49) 


As a conclusion, by (43), (44), (45) and (49), one may sum up the asymptotic behavior of the predictable 
variation process ((A 4 )„)n>o by the following formula. 


(M)n 

y^l-\-a.d 


n^Ai 0 0 

0 A2 0 
0 0 As 


(50) 


It should be noted that the coefficients Ai, 1 < i < 3, are positive because we assume that Voo{x) > 0 
and f{x,t) > 0 in the statement of Theorem 3.3. 


B.3.3 Limit theorems for the vector martingale 

Law of large numbers. We propose to apply the law of large numbers for vector martingales (see 
[21, Theorem 4.3.15]) to the process of interest (Atn)n>o- In the sequel, for any n, Tn denotes the trace 
of the matrix (At)n, while stands for its minimum eigenvalue. First, in light of (50), the trace (7)i)„>o 
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almost surely tends to infinity. Thus we are able to apply the third item of Theorem 4.3.15 of [21] with 
any function a{t) = rj > 0, and we obtain, 





By (50) again, we have, when n goes to infinity, 

Tn = O and 


£n ~ min(A 2 , A 3 ) 


As a consequence, using that - 

of large numbers, as n tends to infinity. 


o(n^) whenever ad + /3 < 1 together with (50), we have the law 



n 


(51) 


Central limit theorem. We now investigate the asymptotic normality of the vector martingale 
{Mn)n>o- We apply Corollary 2.1.10 of [21] with the sequence (a„)n>o defined by a„ = The 

first assumption of this result is obviously satisfied: the sequence (a^^(A4)n)n>o almost surely converges 
to some positive semi-definite matrix. Indeed, by (50), 


-l — ad—i3 


(Mk 


l+ad+p " ^ 

0 0 0 

0 0 0 


E{x,t,a,l3), 


(52) 


where E(a;, t, a, P) is a degenerate variance-covariance matrix. As a consequence, we only have to check 
Lindeberg’s condition in order to establish the central limit theorem for {Mn)n>o- In other words, we 
have to prove that, for any e > 0 , 


1 


a l+Cl(i+/3 


Ee 


i=i 




1+ad+P 

|Afj-Mj_i|>en 2 


F,-l 


0 . 


(53) 


For any 1 < j < n and 1 < fc < 3, we study the three components We have 






U + 1 ) 


nd+iS 


v^wo 


Zj - X 


Ki 


Sj+i — t 


■I 


IR^ xIR 

ad 


Q{^{Zj-i,Sj),x + uvj)f{x -I- uvj,t + vwj) IKd(u)IKi(u) Ad(du) du. 


U + 1 ) 




- / Q{^{Zj-i,Sj),x + uvj)G{x + uvj,t)Kd{u) Xd{du), 

jRd 

[ Q{^{Zj-^,Sj),x + uvj)Kd{u)\d{du). 
^0 V Vj J j^d 


3 3 

Thus, we obtain 




{n + l)°“^+^ 


VglUO 


^d||oo||IKi||oo + ||Q||oo||/||oo, 


vg 


/ l+ad+0\ 

Together with the condition ad + < 1, \Mj — Mj-i\ = oin 2 1 . As a consequence, there exists 

an integer such that for any n > the event < \Mj — Mj-i\ > en 2 Ms almost surely empty. 
This shows Lindeberg’s condition (53). Finally, by (52) and (53), we obtain, as n tends to infinity, 

l+ad+S ^ 

n 2 ^ A/* (O 3 , E(a:,t, a,/?)). (54) 
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C Proof of Proposition 3.14 


By virtue of the ergodic theorem (see Proposition A.2) applied to the Markov chain {Zn, Sn+i)n>o, we 
have, as n tends to infinity, 


>n,n,p / \ CL-s 


(a) 




TT 2 


/ 


e"(C,^'x(C))l(e,(5),+oo)(*)Moo(dC X dt), 


conditionally to o-{Zo, Si, ..., Z„-i, S„). Together with the definition of fj,ao (31), the expression (9) of 
S{x, •) and Remark 2.2, we obtain 

f (55) 

JC:c TT 2 pd-l Jj,c,p 

conditionally to cr(Zo, ^i, ..., Zn-i, Sn). By definition (24) of lx,p we have 

[ g"(c.04e))G(e,e.(c))^oo(c)A,(dc)= / If g”(?,0,(c))G(c,0,(c))*^oo(e)dc)Arf_i(d2/). (56) 


In addition, one has 


r* (y) 


f e"(C,e.(C))G(C, 0 ^(O)^oc(C)de= f r($(y,-t),t)G(-I>( 2 /,-t),t)i 2 ^(-I>(y,-t))|Vt$(y,-t)|dt. 

JCy Jo 


As a consequence, under Assumptions 3.13, a conscientious calculus together with eq. (25) shows that 


L 


L 


G^{^,ey{0)G{^,ey{0>^i0d^- / e"(^,r4C))G(C,r40)j2^(e)dC 


= 0{\x-y\). 


(57) 


We obtain the expected result from (55), (56) and (57) and by remarking that the normalizing constant 
is Ad—i(Dx,p). 

Remark C.l In order to prove (57), one may split the integral of interest into two terms: an integral 
between 0 and t~ (x) At~ (y) and a remainder term integrated between t~ (x) At~ {y) andt~ {x)yt~ (y). The 
first one is clearly upper bounded by the integral of a Lipschitz function between 0 and t~ {x) + \x—y \ [t~]Lip, 
thus bounded by an integral between 0 and t~ (x) + p[t~]Lip- The second integral is obviously bounded by 

\\Q II oo II r^oo II oo 11 II oo [t ]i/ip|^ y\• 
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Figure 2: Two representations of the same simulated path of the TCP-like model of interest until the 
jump. A vector field graph is given on the left, while we observe the trajectory of the hrst component versus 
time on the right. 



Figure 3: Estimation of the invariant distribution of the post-jump locations computed from the 20 000 first 
jumps of the TCP-like model. 
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Figure 4: The choice of the bandwidth parameter a appearing in the criterion /?" is obtained by minimizing 
the cross-validation estimate of the related ISE. The parameter p seems to have only a small 

influence over the minimization: the estimated error is computed from p = 0.005, p = 0.01 and p = 0.02 
from left to right. 
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Figure 5: The optimal index Q is calculated by maximizing the criterion ^ € Cx, computed with 

the optimal parameter a. On the left side of the figure, we compare the estimate k”(C) (full line) with 
F^(5)G'(5,T 2 ,(^)) (dashed line) both computed from a = 0.1, which confirms that the estimation performs 
pretty well. This quantity seems to admit one and only one absolute maximum. 
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Figure 6: The parameter a plays a crucial role in the estimation of We compare here the curves 
computed from a = 0.1 and already presented in Figure 5 (full and dashed lines) with the too oscillating 
estimate obtained from a = 0.4 (dotted line). 
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Figure 7: The optimal bandwidth parameters a and (3 implicitly appearing in 1F”(C) ^^(C))) ^ ^ Cx, are 
obtained by minimizing the cross-validation estimate (a,/3), computed from pi = p 2 = 0.1 here, of 

the related ISE. The parameter /3 seems to have only little influence over the estimation error in comparison 
to a. 
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Figure 8: The estimator X^{x) is computed from different values of the index ^ € Cx, and from the optimal 
bandwidth parameters and (a'^,/3'^). Over the 100 replicates, the optimal points are most of the 
time located at 0.55 and more than 90 times between 0.5 and 0.6 (enhanced by gray colors), which seems to 
correspond with the estimators with least bias and variance. In particular, we obtain a better result than 
by maximizing the estimated invariant measure (index around 0.35), see Remark 3.11. 



Figure 9: Simulated trajectory of the bacteria in the unit disk until the 100*^ jump. 
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Figure 10: For each of the 9 target points {xi,X 2 ), the boxplot of the optimal estimators A2 ,Axg) for 9 
taken in a uniform discretized grid of [0, 27r] with step tt/S is presented. The black thick lines in the boxplots 
correspond to the aggregated estimate A(a:i,a; 2 ) defined in (30), while the dashed line represents the true 
jump rate. 
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Figure 11: For each of the 9 target points, the curve ^xi,x -2 as well as the estimate \{x\,X 2 ) have been 
computed from n = 100 000 observed data and are provided (left). The boxplot of these 9 estimates is also 
given (right) and presents no bias and a small dispersion. 
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Figure 12: For each of the 9 target points, the curve as well as the estimate A(xi,a; 2 ) have been 

computed from n = 20 000 (top) and n = 50 000 (bottom) observed data and are provided (left). The 
boxplots of these 9 estimates are also given in both cases (right). 



Figure 13: For each of the 9 target point {xi,X 2 ) indexed from 1 to 9, estimated jump rates A(a;i,a; 2 ) 
computed from different datasets. 
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Figure 14: Boxplots of the 9 estimated jump rates A(xi,a: 2 ) computed from datasets of different sizes. 
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Figure 15: Virkler’s dataset contains 68 independent crack growth histories starting from qq = 9 mm. 
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Figure 16: Material parameters m and log(C) are strongly linked by a linear relationship used to reduce the 
dimension of the underlying model: log(C') = —9.25 — 5.89 x m + e. 
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Parameter m Crack length 

Figure 17: Estimation of the criterion Ka{m) for different values of the target length a: a = 25mm (top, 
left), a = 35mm (top, right) and a = 45 mm (bottom, left). The relationship between the optimal parameter 
m maximizing Ka{m) and the target crack length a is also presented (bottom, right). 
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Figure 18: Estimation of the jump rate A(a) for different crack lengths a between 25 mm and 45 mm in the 
stochastic model of fatigue crack propagation from Virkler’s dataset. 
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